Landau-Zener sweeps and sudden quenches in coupled Bose-Hubbard chains 



C. Kasztelan\ S. Trotzky^'^ Y.-A. Chen^'^, I. Bloch^-^, I. P. McCuUoch*, U. Schollwocki^^ ^^^^ q Qtso^^^ 
^ Physics Department, Arnold Sommerfeld Center for Theoretical Physics and Center 
for NanoScience, Ludwig-Maximilians-Umversitat Miinchen, 80333 Miinchen, Germany 
^ Physics Department, Quantum Optics Chair, Ludwig-Maximilians-Universitdt Miinchen, 80799 Miinchen, Germany 
^ Max-Planck Institut fiir Quantenoptik, 85748 Garching, Germany 
* School of Physical Sciences, The University of Queensland, Brisbane, QLD 4 072, Australia 
^ Institute of Advanced Study Berlin, 14193 Berlin, Germany and 
^ Laboratoire Matenaux et Phenomenes Quantiques, Universite Paris 
Diderot-Paris 7 and CNRS, UMR 7162, 75205 Pans Cedex 13, France 
(Dated: January 13, 2013) 

We simulate numerically the dynamics ol strongly correlated bosons in a two-leg ladder subject 
to a time-dependent energy bias between the two chains. When all atoms are initially in the leg 
with higher energy, we find a drastic reduction of the inter-chain particle transfer for slow linear 
sweeps, in quantitative agreement with recent experiments. This effect is preceded by a rapid 
broadening of the quasi-momentum distribution of atoms, signaling the presence of a bath of low- 
energy excitations in the chains. We further investigate the scenario of quantum quenches to fixed 
values of the energy bias. We find that for large enough density the momentum distribution relaxes 
to that of an equilibrium thermal state with the same energy. 

PACS numbers: 37.10. Jk, 05.70.Ln 



The investigation of dynamical phenomena in quantum 
many-body systems provides an efficient way to probe 
the physical properties of the underlying ensemble and 
to tackle fundamental questions of statistical mechanics. 
Of special interest in this regard is how the quantum 
dynamics of two coupled modes is changed in a many- 
body setting [IHSl El] and how a closed quantum system 
far from equilibrium does or does not equilibrate (see 
Refs. [5H5] and references therein). The high degree of 
control achieved in current experiments with ultracold 
atoms in optical lattices [3] allows for clean and quanti- 
tative studies of quantum many-body dynamics jlOH16j . 
Recently, a generalization of the famous Landau-Zener 
(LZ) sweep [T7] to a pair of coupled one-dimensional 
quantum gases of strongly correlated bosonic particles 
was addressed experimentally [18] . One intriguing result 
of this experiment was the observation of a breakdown of 
adiabatic inter-chain transfer in slow sweeps of the inter- 
chain bias. 

In this Letter, we present a numerically exact study 
of the LZ dynamics in two coupled Bose-Hubbard chains 
subject to a time-dependent inter-chain bias energy, us- 
ing the time-dependent density matrix renormalization 
group (i-DMRG) method [H]. We show that the break- 
down of adiabatic transfer as observed in Ref. [IH] in slow 
sweeps far away from the ground-state is always accom- 
panied by a dramatic broadening of the quasi-momentum 
distribution of atoms in the two legs. This provides 
strong evidence that the responsible mechanism for the 
breakdown is the coupling to an internal bath of low- 
energy momentum excitations. Finally, we study the 
quantum dynamics emerging after sudden quenches of 
the bias energy where the same mechanism causes a re- 
laxation to a steady state similar to an equilibrium ther- 



mal state with the same energy. 

Setup and model. We consider a two-leg ladder formed 
by two coupled Bose-Hubbard chains [see Fig. [ija)] de- 
scribed by the Hamiltonian 
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where fo^.o- (i — \...La, a = L, R) are the local an- 
nihilation operators and ni^„ — S^^Si o-. The parame- 
ters J|| and J are the intra-chain and inter-chain tunnel 
couplings, respectively, \J is the on-site interaction en- 
ergy, and A(f) is a time-dependent energy bias between 
the two chains. In the following we fix J = /i = 1 
and we call n = N/Lg the density of the system with 
N — ^{fii^a)- The zero-temperature properties of the 
model ([I]) have been investigated recently using static 
DMRG and mean-field methods [20, . 

We initialize the system in the ground state of Hamil- 
tonian Eq. ([T]) with a large positive bias A(t = 0) — )■ -l-oo 
so that all particles are initially located in the left leg. 
This ground state is found from a static DMRG calcula- 
tion. We distinguish two different scenarios for the time- 
dependence of the bias A{t): 

A(i>0) = |^° + "'' (linear- «^eep) 

1 A/ , (sudden quench) . 

The first case corresponds to a linear LZ sweep Aq — > 
— Aq of the bias with constant rate a — — 2Ao/T, T 
being the sweep time. The second case amounts to a 
sudden change (quench) of the bias to a constant value 
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FIG. 1. (Color online) (a) Sketch of the two-leg ladder de- 
scribed by the Hamiltonian Eq. (1) (b) Instantaneous right- 
leg population n^{t) during ground-state (dashed lines) and 
inverse sweeps (solid lines) with various rates a, calculated 
with n = 0.75 and La = 8. The Hamiltonian parameters are 
|Ao| = 18.2, J|| = 0.38, and U = 1.58. (c) Transfer efficiency 
nR,(r, 2TT/a) for inverse sweeps with the same parameters and 
n = 0.75, 1, 1.25 (solid lines). The dashed line is the two-level 
LZ probability phz{a)- Open circles represent experimental 
data taken from Ref. [18] for the corresponding set of control 
parameters [see Appendix B]. The inset shows the maximum 
transfer efficiency as a function of the density for Ls = 8 
(open circles) and 12 (filled circles) and the linear extrapola- 
tion to higher densities (dashed and solid line, respectively). 
The blue line marks the maximum reached in the experiment. 

Af. We use i-DMRG [H] to calculate the exact quantum 
dynamics in both situations for system sizes up to Ls ~ 
16. In all the simulations we fix the truncation error to 
e — 10~^ and check that all results have converged. 

Linear sweeps. We first discuss LZ sweeps of the en- 
ergy bias, focusing on inverse sweeps, where the left leg 
is initially higher in potential energy than the right leg 
(Ao < 0, a > 0). The corresponding counterpart we 
call a ground-state sweep (Ag > 0, a < 0). In the limit 
Ls = 1 and U — = 0, the problem reduces to the 
original two-level LZ problem |17j . where the probabil- 
ity to reach the right-leg as a function of the sweep rate 
Plz(q;) = 1 — exp(— 27r/a) is the same for both. 

We consider the fraction riR(t) = of 
particles in the right leg at a given time t. This is plot- 
ted in Fig. [T|b) for both types of sweeps and three ab- 
solute values of the sweep rate a. In Fig. [TJc), we plot 
the final transfer efficiency n^{T) reached at the end of 
the inverse sweep as obtained from a number of such 
traces with varying sweep rates [see Appendix A]. We 
compare the results for n = 0.75, 1, 1.25 to experimental 
data taken from Ref. [15] for the respective set of Hamil- 



tonian parameters [see Appendix B]. In the experiments, 
the effective density was significantly higher as will be 
discussed in detail below. 

For fast sweeps, corresponding to |a|/27r ^ U, Jy, the 
curves for the ground-state and for the inverse sweep lie 
on top of each other [Fig. [TJb)-z]. This is due to the fact 
that practically no intra-chain processes can occur during 
the sweep and hence the problem is effectively reduced 
to the original two-level LZ problem. As a consequence 
the DMRG results in Fig. [ijc) collapse onto plz (a) for 
small values of 27r/a. 

Differences between the two kinds of sweeps become 
evident as a decreases and intra-chain processes start 
to affect the dynamics. Here, we find the transfer effi- 
ciency for the inverse sweeps to be reduced relative to 
the ground-state sweep [Fig. [lJb)-H]. This can in parts 
be understood in the limit of isolated double wells filled 
with interacting particles (J|| ^ 0, U > 0). Here, the 
dynamical problem can easily be solved numerically by 
direct integration of the Schrodinger equation and even 
analytically when assuming U ^ 1 [TH • The analy- 
sis of this scenario shows that with increasing number of 
particles the transfer efficiency is enhanced in the ground- 
state sweep, while it is reduced in the inverse sweep [see 
Appendix C]. 

For slow enough sweeps, corresponding to |a|/27r < 
min(J7, J||), the transfer efficiency for the inverse sweep 
starts to decrease as the rate a decreases, as it was found 
in Ref. [THj. This is in sharp contrast with the behavior 
for the ground-state sweep, where the transfer efficiency 
approaches unity [Fig. [T|b)-Mi]. We emphasize that the 
breakdown of adiabaticity in the inverse sweep cannot be 
obtained within the isolated double- well picture and is di- 
rectly related to the possibility for the system to generate 
excitations along the chains. Since these excitations have 
to be created via collisions, the breakdown is more effec- 
tive in systems with higher densities, as it is evident in 
Fig. [ij^c). Eventually, for infinitely slow sweeps {a — >■ 0) 
we expect the transfer efficiency to rise to 1 again due to 
the finiteness of the system. 

We note, that the experimental results taken for com- 
parison were obtained from a two-dimensional array 
of pairwise coupled chains with inhomogeneous density 
[18) . In each of these ladders the maximum fidelity 
is reached for different sweep rates and takes different 
values. Therefore, the maximum in n^(T) is less pro- 
nounced than in the numerical results for a single ladder. 
From ground-state DMRG calculations respecting the ex- 
perimental geometry, we find that the average density in 
the center of the chains was {n{z = 0))av — 2, whereas 
the overall average density was (n)av — 1-4. A rapid 
growth in entanglement entropy [21] |22] prevents us to 
access the dynamics with these densities - or of the full 
inhomogeneous system - directly in the simulations. In 
the inset of Fig. [ijc) we plot the maximum efficiency 
JT-R.max reached in the DMRG simulations as a function 
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FIG. 2. (Color online) Broadening of the left- and right-leg 
momentum distributions during the inverse sweep: (a) Quasi- 
momentum distributions Uka of the left and the right leg as 
a function of the reduced sweep time 27r/a as obtained from 
DMRG for n = 0.75 and La — 16. The arrows mark peaks 
in the momentum distributions at = ±7r which appear for 
27r/a > 1. (b) The corresponding quasi- momentum width 
{fc^)a- of both legs (filled symbols) plotted together with the 
results for Ls — 8 (open symbols) where slower sweeps are 
accessible. The dashed vertical line marks the sweep rate for 
which riR, reaches its maximum [see Fig. Wic)] . The Hamilto- 
nian parameters are the same as for Fig.TTl 



of n and for Ls ~ 8, 12. We find riR.max to depend ap- 
proximately linearly on n and only weakly on L^- Lin- 
ear extrapolation of our numerical results for nj^ to 
higher densities yields a crossing with the experimentally 
measured value at n ~ 1.8. 

Having recovered the downturn of the transfer efh- 
ciency observed in the experiments, we now turn to 
the discussion of the experimentally accessible quasi- 
momentum distribution of atoms in both legs |9I ■ The 



latter is defined as rn^a — L ^ J2r. 
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where the sum extends to all lattice sites and k is given 
in units of the reciprocal lattice constant. In Fig. [2ja), 
we show Tifco- for the left and the right leg as a function 
of 27r/a and k and their respective widths in Fig. [2|^b). 

For very fast sweeps, the left- and right-leg distribu- 
tions are equal up to a constant factor and therefore 
have equal widths. Since intra-chain processes are ab- 
sent on this timescale, each individual Bloch state un- 
dergoes a separate two-level LZ transition. Therefore, 
the momentum distribution in the right leg reduces to 
nkK = P'L7.{ct)nkij{t = 0). As the sweep rate decreases, 
the distributions of both legs show a completely different 
behavior. In the left leg, the peak around k — Q \s more 
and more depleted, whereas new local maxima emerge 
near k — ±7r, signaling that the motion of atoms on 
neighboring sites in the chains becomes correlated. The 
width of the quasi-momentum distribution grows 

correspondingly, reaching a maximum value shortly be- 
fore the transfer efficiency becomes maximal (dashed ver- 
tical line). 



For even slower sweeps, the peaks at k = ±7r in the 
left leg spread out and make the quasi-momentum dis- 
tribution appear more uniform. This, together with a 
counterflow of low-momentum components from the right 
to the left leg, causes the width (fc^)L to decrease. The 
quasi-momentum distribution of the right leg, in con- 
trast, shows the condensate peak at fc = but the peaks 
at the edge of the Brillouin zone are almost absent. The 
associated width {k"^)^, increases about linearly for the 
full range of values 27r/a covered by our simulations. 
These results provide numerical evidence for the coupling 
to low-energy momentum modes being the responsible 
mechanism for the breakdown of adiabaticity in the in- 
verse sweeps. 

Quantum quenches. A more detailed "stroboscopic" 
understanding of this decay mechanism can be obtained 
from the quantum quenches denoted by the second case 
m Eq. ^. Since the Hamiltonian is time- independent for 
i > 0, the energy E of the system is a constant of motion 
throughout the subsequent time-evolution. In Fig. [sj^a) 
and (b), we plot the instantaneous values of the fraction 
nii{t) of particles and the width of the right-leg quasi- 
momentum distribution (A;^)R(t), respectively, for three 
different values of the bias A f where the system is again 
initialized with large and positive. Both observables 
show oscillations which dampen out towards a non-zero 
value which depends of A/ . 

In Fig.jsjjc), we plot the long-time value of the right-leg 
population = ^^(t — > c») as obtained from the time 
traces for different values of A/. The dots and the error 
bars represent, respectively, the average value and the 
amplitude of the last oscillation that we simulate numer- 
ically [see Appendix D] . As the bias A / crosses zero and 
becomes negative, reaches a maximum value slightly 
above 0.5 and then decreases. Any loss of potential en- 
ergy due to the transfer of atoms to the empty chain 
must be counterbalanced by an equal increase in kinetic 
energy. Since this condition can hardly be satisfied when 
the detuning is large and negative, more and more atoms 
remain self-trapped in the left well (shaded region). 

For comparison, we performed static DMRG simula- 
tions for an interacting Bose gas on the two-leg ladder 
with bias Ay at thermal equilibrium. We chose the tem- 
perature in these simulations such that the total energy 
matches with the conserved energy E in the quench. The 
respective results for [solid line in Fig. [sj^c)] are in 
very good agreement with the long-time values of the 
dynamical evolution. On the other hand, the same cal- 
culation for an ideal Bose gas with E = —2J\\nLg (dashed 
line) shows significant discrepancies, indicating the cru- 
cial importance of interactions. The inverse (effective) 
temperatures /3 found from the matching condition are 
plotted as inset in Fig. |3|c). At the point A/ = — 4J||, 
/? becomes zero irrespective of the density or the interac- 
tion strength. This implies that quantum quenches with 
Ay < — 4J|| cannot be associated to an equilibrium ther- 
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FIG. 3. (Color online) Quantum dynamics after a sudden 
quench: (a) Time-traces of the right-leg density nnit) and (b) 
the width {k^)ii of the momentum distribution in the right leg 
for n — 1.5 and different values of A/, (c) Long-time value 
riR of the right-leg density versus Af (n = 1.5). Points and 
error bars refer to the time-averaged value and amplitude of 
the last oscillation, respectively [see Appendix D]. The solid 
line denotes the right-leg population as obtained for an in- 
teracting thermal Bose gas with the same total energy as the 
initial state in the quenches. The dashed line is the prediction 
for an ideal Bose gas. The shaded area indicates A/ < — 4J||. 
The inset shows the inverse temperatures used in the finite- 
temperature calculations, (d) Long-time value of the width of 
the momentum distribution versus detuning for various den- 
sities (circles). The solid and dashed lines are the results of 
the finite temperature calculations as in (b). For all plots, the 
size of the system is Lg — 16 and the Hamiltonian parameters 
are as in Fig. [l] 



mal state with a finite positive temperature. 

In Fig. [Sf^d) we plot the long-time value of the width 
{k^)'R of the right-leg quasi-momentum distribution as a 
function of A^: and for different densities. Within our 
numerical accuracy, the quasi-momentum widths in the 
left and in the right legs coincide for long times, unless 
Af < — 4J||, where the dynamics are found to be ex- 
tremely slow. We again compare the long-time values 
with the results of the finite-temperature DMRG calcu- 
lation and find good agreement for the largest density 
(n = 1.5), especially away from Af = 0. For lower densi- 
ties, and closer to A/ = 0, we observe a less pronounced 
broadening of the momentum distribution. Here, colli- 
sions are much less effective and a complete equilibration 
cannot be reached before finite-size effects become im- 
portant [see Appendix D]. At Af = — 4J||, where the 
effective temperature diverges, the momentum distribu- 
tion is flat, corresponding to (/c^) = 7r^/3 for any finite 
density. 

Together with the equilibration of right-leg density, our 
findings suggest that the system relaxes towards an equi- 
librium state close to the Gibbs-ensemble associated with 
the initial state's particle number and total energy. The 
responsible decay mechanism is the same that leads to the 



breakdown of adiabaticity and the accompanying broad- 
ening of the momentum distribution for slow inverse LZ 
sweeps. 

In conclusion, we have studied numerically both LZ 
sweeps and quenches in a two-leg ladder system of 
strongly interacting bosons. For inverse sweeps, we have 
recovered the breakdown of adiabatic transfer of parti- 
cles to the initially unoccupied right leg as it was re- 
cently found in experiments |I8j. We have shown that 
this phenomenon is preceded and accompanied by a fast 
broadening of the momentum distribution of the initially 
filled left leg, providing numerical evidence for the cou- 
pling to an inner bath of low energy momentum excita- 
tions. Finally, we have investigated the underlying decay 
mechanism by studying quantum quenches of the energy 
bias. We have found strong evidence for the system to 
approach a thermal state whose temperature is set by 
the bias. Our findings provide detailed insight into the 
dynamics of strongly correlated quantum many-body sys- 
tems in low dimensions which can be probed in current 
experiments with ultracold atoms. 

This work was funded by the DFG (FOR 635, FOR 
801) and the EU (NAMEQUAM). 
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APPENDIX A - TIME-DEPENDENT DMRG 
SIMULATIONS 



From our simulations we seek to obtain the final trans- 
fer efficiency riR in the limit of |Ao|,T — cx) in order 
to compare or findings to the original LZ problem [17]. 
However, for finite values of |Ao| and T, when fixing the 
bias to — Aq for t > T, the right-leg population is not con- 
stant, but oscillates around an average value. This fact 
complicates the extraction of a final value for the transfer 
efficiency especially for fast sweeps and is caused by two 
reasons. On the one hand, the initial state calculated at 
A{t = 0) = 100 with all particles located on the left leg 
is projected onto the eigenstates of the Hamiltonian (1) 
with A = Ao (|Ao| < A{t = 0)). In other words, the 
initialization of the LZ sweep is already practically non- 
adiabatic. On the other hand, for finite sweep rates, the 
sweeps are not fully adiabatic. The final state hence is a 
partial projection of the initial state onto the eigenstates 
of (1) with A = — Aq. Since these eigenstates consist of 
contributions from both legs for values of |Ao| not much 
larger than 1, coherent population oscillations will occur 
for t > T. 

We therefore rescale both the energy bias Aq — tAq 
and the sweep time T' = rT by the same factor r, which 
leaves a = 2Aq/T unchanged but makes the simulation 
more adiabatic. In practice in our numerics we choose 
this 'multiplication' factor between r — 7 for fast sweeps 
and r = 2 for slow sweeps. For each final value of tir in 
Fig 1(c), we take the average over the, now, much smaller 
residual oscillations at the end of every single sweep. 



APPENDIX B - CHOICE OF HAMILTONIAN 
PARAMETERS 

In the experiments of Ref. [18J, a superlattice of 
the form Vx{x) = VxsSm'^{2Trx/Xa:s) + Vxism'^{2Trx/Xxs) 
was created along the x-direction by superimposing two 
standing- wave laser fields with wavelengths X^s = 765 nm 
and Xxi = 1530 nm. Additional transverse lattices with 
wavelengths Xy^^ = 844 nm and respective depths Vy^z 
completed the setup. For all simulations presented in 
this Letter, we fix the lattice depths to be Vxs — 15 -Efoci 
Vooi = 40 = 30 Ef^^ and 14 = 4 E^^^ in accordance 

with a specific measurement in Ref. [18] and calculate 
the corresponding Hamiltonian parameters from the sin- 
gle particle bandstructure. All lattice depths are given in 
units of their respective recoil energies E^^^ = h"^ / {2mXi). 




FIG. SI. (Color online) Transfer efficiency for a ground-state 
sweep in a double-dot. (a) Plot of the analytical result for tir 
in the Hmit of (7 ^ 10 for n = 1,2,3 bosons in the double 



well in the case of a ground-state sweep (solid lines, Eq. (SI l) 
and an inverse sweep (dashed lines, Eq. (S3l). (b) Plot of 



the transfer efficiency as obtained from direct numerical inte- 
gration of the Schrodinger equation for [/ = 1 and the same 
conditions as in (a). 



APPENDIX C - LZ SWEEPS WITH FEW 
PARTICLES IN A DOUBLE WELL 

As long as Q;/27r > J|| in the LZ sweeps, the dynamics 
of distant rungs in the ladder is still uncorrelated. In 
this limit, it is helpful to consider the transfer efficiency 
for a LZ sweep within an isolated double well occupied 
by n interacting particles, formally corresponding to a 
"ladder" with = 1- 

Ground-state sweeps. For a strong repulsive on-site 
interaction U ^ n, the ground state undergoes TV in- 
dependent single-particle tunnel resonances, as the bias 
A is changed [231 [21] ■ At each of these resonances, two 
Fock-states \n — j, j) and jn — j — 1, j -I- 1) {j — Q . . .n—1) 
are coupled and the corresponding coupling matrix ele- 
ment for the bosons is —Jn,j = + 1) < — 1. 
The propbability to diabatically cross one of these reso- 
nances is given by pnj{a) = exp(— 27r(n — + l)/a). 
If the j-th resonance is crossed diabatically, the double- 
dot remains in the state \n — When calculating 
the transfer efficiency of a sweep, we can neglect any 
higher order crossings between the excited states pos- 
sibly reached, since the associated coupling matrix el- 
ements are typically suppressed as U~'^ , where v is the 
order of the respective resonance. Thus, the final transfer 
efficiency eventually becomes 



1 - 



(SI) 



Inverse sweeps. Other than the ground state in the 
double well, the highest excited state undergoes only a 
single tunneling resonance |n, 0) — >■ |0, n) at A = 0, where 
all particles have to change sides simultaneously. In the 
limit of strong repulsion, this tunnel process is suppressed 
as J/^"" ^ 1 (rt > 1). Therefore, much slower sweeps 
than in the ground state are needed to achieve full trans- 
fer. If this resonance is crossed diabatically and the sys- 
tem remains in state |n, 0) , it will arrive at another tunnel 
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resonance to |1, n — 1) which is of (n — l)-th order, as A is 
increased. Further along this path, the system will face 
a total of n tunnel resonances of decreasing order and 
thus for J7 3> 1 increasing coupling. The coupling matrix 
element for the v-th order resonance |n, 0) — ?> |n — v, v) 



in leading order is 

The probability to follow the eigenstate adiabatically at 

resonance v is thus 1 — 'p'r? {p) = 1 — exp[— 27r J,!'^'' /a]. 
Since this adiabatic crossing results in a final right-well 
population the overall transfer efficiency becomes 



n — 1 



nR(a) = (1 + — - + ... + - (a) . . .p(f)(a)(l - p(i)(a)) 



^ n— 1 n 



(S3) 



fj,—0 f—n—^ 



In Fig. SI |a), we plot the transfer efficiency for both a 



ground-state sweep and an inverse sweep with n = 1, 2, 3 
particles and U — 10. In the case of the ground-state 
sweeps, the enhancement due to bosonic statistics at each 
resonance leads to an increase of the transfer efficiency 
with n. For the inverse sweeps, however, the transfer effi- 
ciency is reduced with increasing particle number due to 
the higher order effective coupling matrix element. Only 
the single-particle transfer |n,0) — > |rt — 1, 1) is of first 
order, so that each of the curves nearly follows the single- 
particle result until a value of from where on the 
transfer efficiency is significantly reduced. This qual- 
itative behavior survives for much weaker interactions 
from a direct numerical integration of the time-dependent 
Schrodinger equation [see Fig. |Sl[b)]. 



APPENDIX D - LONG-TIME VALUES FROM 
THE QUENCHES AND FINITE-SIZE EFFECTS 

To extract the long-time values from the time-traces 
for nYi{t) and {k"^) in the quench scenario [Fig. 3(c) and 
(d)], we average the traces over the last oscillation pe- 
riod accessible in our t-DMRG simulations. For a den- 
sity of n = 1.5, these residual oscillations have only 
small amplitude and no significant drift is observed any- 
more [see Fig. 3(a) and (b)]. For low densities, i. e. 
n = 0.75, we are able to follow the dynamics up to times 
t > T* = L^/(2J||) ~ 21 {Ls = 16, J|| = 0.38), where 
finite-size effects caused by reflections at the boundaries 
of the system become important. Here, we estimate the 
"asymptotic value" by taking the average over the oscil- 
lation period located around i = 15 < T*. Since colli- 
sions are rare for such low density, the oscillations am- 
plitudes are slowly damped and the system cannot relax 
fully within the time T* , leading to the pronounced de- 
viations from the results for the thermal state shown in 
Fig. 3(d). 



